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There are inherent uncertainties and errors associated with using Computational Fluid 
Dynamics (CFD) to predict the flow Held and there is no standard method for evaluating 
uncertainty in the CFD community. This paper describes an approach to validate the 
uncertainty in using CFD. The method will use the state of the art uncertainty analysis 
applied to the ke-realizable turbulence model to predict the velocity uncertainty of a 
backward facing step. 


Nomenclature 


£21 = solution changes medium to fine grid 

£32 = solution changes coarse to medium grid 

e a 2 1 = extrapol ated error 

GCI fi ne 21 = grid convergence index 
h = representative grid size 

p = observed order 

Rk = convergence parameter 

f* 2 i = ratio of grid sizes between grid 1 and 2 

1*32 = ratio of grid sizes between grid 3 and 2 

Ski = solution variable for fine grid 
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Sk2 = 

Sk3 


s 21 = 

3 ext 


Sl = 

Su = 

Uosdllatory 

Umonotonic 


solution variable for medium grid 
solution variable for coarse grid 
extrapolated solution variable 
lowest solution variable 
highest solution variable 

= uncertainty for oscillatory portion of the solution 
= uncertainty for monotonic portion of the solution 


I. Introduction 

is the current state of the art and industry standard used for flow field predictions and analysis; however 
CFD has many challenges. There are inherent uncertainties and errors associated with using CFD to predict the flow 
field, and there is no standard method for evaluating uncertainty in the CFD community l . 

Some potentials errors include physical approximation error, computer round-off error, iterative convergence 
error, discretization errors, computer programming errors, and usage errors 2 . An uncertainty, as defined by the 
American Institute of Aeronautics and Astronautics (AIAA), is a potential deficiency in any phase or activity of 
modeling and simulation that is due to the lack of knowledge 3 . An example of an uncertainty in performing a CFD 
analysis is turbulence modeling 4 . There is a lot about turbulence modeling that is not understood 4 . There has been 
progress in estimating the uncertainty of CFD, but the approaches have not converged V 

A thorough literature review has been performed to determine the best method to evaluate the uncertainty in 
CFD predictions. Both major journals in mechanical and aerospace engineering, AIAA and ASME, have published 
articles on this subject. The ASME method has been adopted by many researchers and provides a detailed approach 
to calculate uncertainty in CFD from different levels of grid refinement. The method published by the ASME 
Journal of Fluids Engineering is the state of the art for determining the uncertainty in CFD predictions and will be 
used for the proposed research problem. 

A CFD model has been created of a backward facing step using ANSYS FLUENT and OpenFOAM. The 
backward facing step was solved previously for re-attachment length by Celik and Karatekin 5 . The backward facing 
step induces turbulence into the flow field and will provide adequate physics to compare the uncertainty of different 


turbulence models using the Comprehensive Approach to Verification and Validation of CFD Simulations 6 . This 
paper provides a detailed uncertainty analysis of the ke-realizable turbulence model for the backward facing step. 

The structure of the paper is as follows. In section 2, the literature review is summarized for CFD uncertainty 
analysis. Section 3 presents the grid refinement study. Section 4 presents the numerical results of the backward 
facing step. Section 5 is the discussion. Section 6 is the conclusion. 

II. Literature Review 

A literature review was performed to determine the “State of the Art” method for calculating CFD uncertainties. 
CFD is extensively used in industry, government, and academia to design, investigate, operate, and improve 
understanding of fluid physics 3 . The rate of growth in using CFD as a research and engineering tool will be directly 
proportional to the level of credibility that the simulation can produce 3 . One needs to evaluate the uncertainty in the 
results of a CFD simulation to postulate a level of credibility. In 1986, The American Society of Mechanical 
Engineers (ASME) Journal of Fluids Engineering published a policy statement stating the need for quantification of 
numerical accuracy V Other journals have issued similar statements 7 . These statements lead to research on the best 
method to determine numerical uncertainty. In 1995, Celik and Zhang published “Calculation of Numerical 
Uncertainty Using Richardson Extrapolation: Application to Some Turbulent Flow Calculations” which used 
Richardson’s Extrapolation method to estimate the uncertainty in CFD 8 . In 1997, Roache published “Quantification 
of Uncertainty in Computational Fluid Dynamics” 1 . Roaches research also used the Richardson Extrapolation 
method to quantify CFD uncertainties. 

In 1998, the AIAA has published a “Guide for the Verification and Validation of Computational Fluid Dynamics 
Simulations” 3 . This document provides guidelines for assessing credibility via verification and validation 3 . The 
document does not recommend standards due to issues not yet resolved, but defines several terms 3 . “Uncertainty is 
defined as a potential deficiency in any phase or activity of the modeling process that is due to lack of knowledge 3 ” 
“Error is defined as a recognizable deficiency in any phase or activity of modeling and simulation that is not due to 
lack of knowledge 3 ” “Prediction is defined as the use of a CFD model to foretell the state of a physical system 
under conditions for which the CFD model has not been validated Uncertainty and error are normally linked to 
accuracy in modeling and simulation 3 . The guide defines four predominate error sources: insufficient spatial 
discretization convergence, insufficient temporal discretization convergence, lack of iterative convergence, and 


computer programming, but does not make claims about the accuracy of predictions \ The guide emphasizes that 
systematically refining the grid size and time step is the most important activity in verification Once the grid has 
been refined such that the discretization error is in the asymptotic region, Richardson’s extrapolation can be used to 
estimate zero-grid spacing 3 . A sensitivity analysis and uncertainty analysis are two methods for determining the 
uncertainty in CFD 3 . The validation test compares a CFD solution to experimental data 3 . The guide has outlined 
the terms and an overall structure to performing validation, but does not offer a quantitative method. 

In 1999, Stem, Wilson, Coleman, and Paterson, E. G., published Iowa Institute of Hydraulic Research (IIHR) 
Report No. 407 titled "Verification and Validation of CFD Simulations" 9 In 2001, the American Society of 
Mechanical Engineers (ASME) Journal of Fluids Engineering published a “Comprehensive Approach to 
Verification and Validation of CFD Simulations” in an attempt to provide a comprehensive framework for overall 
procedures and methodology 6 . Two papers were published on the subject in Parts 1 6 and Parts II 10 and used the 
methodology documented in IIHR Report 407. Numerical errors and uncertainties in CFD can be estimated using 
iterative and parameter convergence studies 6 . The method uses three convergence conditions as possible in 
estimating uncertainties; (1) monotonic convergence which uses Richardson’s extrapolation, (2) oscillatory 
convergence which uses the upper and lower bounds to estimate uncertainty, (3) divergence in which errors and 
uncertainties cannot be estimated 6 . The literature provides an approach for estimating errors and uncertainties in 
CFD simulations for each of the three cases 9 ’ 6 ’ 10 . The approach uses Richardson’s extrapolation, which is not 
new, however; the method has been extended to use input parameters and correction factors to estimate errors and 
uncertainties 9 ’ 6 ’ 10 . The method examines two sources for error and uncertainty: modeling and simulation. 
Examples of modeling errors include geometry, mathematical equations, boundary conditions, turbulence models, 
etc. Ivn . Examples of numerical errors include discretization, artificial dissipations, incomplete iterative and grid 
convergence, lack of conservation of mass, momentum, energy, internal and external boundary non-continuity, 
computer round-off etc. 4 . The method lacks correlations among errors and assumes these are negligible, which may 
be inappropriate for some circumstances 6 . Additionally, the method provides a quantitative approach for 
determining the iterative convergence uncertainty 6 . Iterative Convergence must be evaluated and is typically done 
by monitoring the residuals order of magnitude drop graphically 6 . For oscillatory convergence, the deviation of a 
residual from the mean provides estimates of the iterative convergence 6 . This is based on the range of the 
maximum S v and minimum S L values 6 For convergent iterative convergence, a curve-fit is used 6 . For a mixed 


convergent/oscillatory, iterative convergence is estimated using the amplitude and the maximum and minimum 
values 6 . A method for confirming validation is presented as compared to experimental data 6 . 

In 2008, the International Towing Tank Conference (ITTC) has published “Recommended Procedures and 
Guidelines - Uncertainty Analysis in CFD Verification and Validation Methodology and Procedures” 11 . The ITTC 
guide was largely based off of the methodology and procedures presented in the ASME Journal of Fluids 
Engineering a “Comprehensive Approach to Verification and Validation of CFD Simulations” 11 . Also in 2008, the 
ASME Journal of Fluids Engineering published a "Procedure for Estimating and Reporting of Uncertainty Due to 
Discretization in CFD Applications” 12 . 

In 2011, the National Energy Technology Laboratory (NETL) conference proceedings held a major section 
related to CFD Uncertainty Calculation 13 . Celik presented “Critical Issues with Quantification of Discretization 
Uncertainty in CFD” 13 . The proceedings were based off of the ASME “Comprehensive Approach to Verification 
and Validation of CFD Simulations” 6 . 

Summary of Literature Review: A thorough literature review has been performed to determine the best 
method to evaluate the uncertainty in CFD predictions. Both major journals in mechanical and aerospace 
engineering, AIAA and ASME, have published articles on this subject. The ASME method has been adopted by 
many researchers and provides a detailed approach to calculate uncertainty in CFD from different levels of grid 
refinement. The method published by the ASME Journal of Fluids Engineering is the state of the art for determining 
the uncertainty in CFD predictions and will be used for the proposed research problem. 

III. Grid Refinement Study 

Convergence studies require a minimum of three solutions to evaluate convergence with respect to an input 
parameter 2 . Consider the situation for 3 solutions corresponding to fine S k i, medium S^, and coarse S k3 values for 
the kth input parameter 2 . Solution changes £ for medium-fine and coarse-medium solutions and their ratio R* are 
defined by 2 : 

£21 = Sk2 - Ski 
£32 = Sk3 - Sk2 


Rk = £21 / £32 


(1) 


Three convergence conditions are possible 2 



Monatomic 

Convergent 


Oscillatory- 

Convergence 


(iii) - Cells with 
Divergence Rk>l 


(i) Monotonic convergence. 0< R,. <1 

(ii) Oscillatory convergence: R k < O' 

(iii) Divergence: R k >l 


The quantity of interest for the backward facing setup is velocity magnitude. Three grids were compared, and 
the convergence conditions were determined for every point in the computational domain. This is accomplished 
through interpolation between the medium to coarse grid and the fine to coarse grid. The velocity magnitude from 
the medium and fine grids are interpolated on to the coarse grid. Then the solutions changes, £ 21 , £ 32 , Rk, and 
convergence conditions are calculated for every point in the domain. Figure 1 shows the different convergence 
conditions inside the computational domain for the grid refinement study. 


Figure 1 : Convergence conditions for a Flat plate - Grid refinement 1 



Figure 2: Velocity Magnitude for Flow over Backward facing step - Coarsest Grid (Structured 1,192,000 cells) 


IV. Numerical Results for Backward Facing Step 


The uncertainty associated with the CFD calculation is the compilation of the elemental errors associated with 
each of the numerical, input, and solver errors. This uncertainty can be calculated using a Data Reduction equation 
the form r = r(X ls X 2 ,. . . Xj) is shown in equation 3, below. 


Ucfd - (£i= 1 |( a; J B i ] + 2 £{=1 £i=i+l {(a*.) (ax J [ B l B k\correlated} + Ei=l {( a *.) P i ]) 


Where, 


Bj = the systematic (bias) error associated with variable Xj, 

(Bj B k ) corTe | atcd = the correlated systematic error between variables Xj and X k , 
and Pi = the random error associated with variable Xj. 


For the calculation, the correlated errors and random errors are neglected and the data reduction equation reduces 
to the following, as shown in equation 4. 

«™=(z;..{(S) ! b?)) /! <«> 


A list of variables for the k-e-realizable turbulence model analyzed is listed in Table 1 . 

Table 1: Uncertainty Variables, Xi 


Type of Variable 

Variables XI 

Value 

Bias Error 

Boundary 

Conditions 

epsillon turbulent mixing length dissipation rate inlet (m2/s3) 

05 

0.5 


k turbulent intensity kinetic energy inlet (m2/s2) 

0.05 

0.05 


pressure outlet (Pa) 

101325 

2 % 


velocity inlet (m/s) 

10 

0.5 

Fluid Properties 

kinematic viscosity nu represents air [0-50-100] deg C 

1.79E-06 

[13.6e-06 -> 23.06e- 
06] 

Grid Size 

Method - Uses Oscillatory Uncertainty 

1,192,000 

1,862,500 

3,311,689 


Numerical 

Method - Uses Richardson's Extrapolation (ASME 5 Step Procedure) - Calculated for Velocity at each 

Cell 

Solver 

Open FOAM (SimpleFoam) vs. Fluent 



Turbulence 

Models 

ke-realiable, kwSST, and SpalartAllmaras 




Expanding the data reduction equation for the listed variables as shown in equation (5) in order from top to bottom. 




Each of the variables was analyzed separately for their elemental error sources. The following plots show the each 
variables and their corresponding uncertainty plot as a function of the percent uncertainty in the CFD Velocity 
prediction. The percent uncertainty is calculated by dividing by the local velocity (ie the uncertainty velocity in 
each cell divided by the velocity in each cell). 


The uncertainty for each of the following was calculated as shown in equation 3 for each cell using the following 
method outlined by Stem, Wilson, Coleman, and Paterson 2 . S is the simulated result. For this case it is the upper 
velocity Su and the lower velocity S L . 


U Oscillatory ~ ? (*$ 1 / $l) 


(3) 


epsilion turbulent mixing length dissipation rate inlet (m 2 /s 3 ) 

For a value of 0.5 +/- 0.5 m /s , the uncertainty in the velocity prediction was 0 — 1. 155 percent as shown in Figure 
3. 



Figure 3. Epsilon Turbulent Mixing Length Dissipation Rate Inlet — Velocity Uncertainty Percentage 





k turbulent intensity kinetic energy inlet (m 2 /s 2 l 


For a value of 0.05 +/- 0.05 m 2 /s 2 , the uncertainty in the velocity prediction was 0 - 0.785 percent as shown in 
Figure 4. 



Figure 4: k Turbulent Intensity Kinetic Energy Inlet - Velocity Uncertainty Percentage 


Pressure outlet (Pa) 


For a value of 101325 +/- 2% Pa, the uncertainty in the velocity prediction was 0-20 percent as shown in Figure 5. 



Figure 5: Pressure Outlet - Velocity Uncertainty Percentage 


Velocity Inlet (m/s) 


For a value of 10 +/- 0.5 m/s, the uncertainty in the velocity prediction was 0 - 6.558 percent as shown in Figure 6. 



Figure 6: Velocity Inlet - Velocity Uncertainty Percentage 


iflnematic viscosit y nu=17.06e-06 [13.6e-06 -> 23.06e-06l (m 2 /s) represents a i r fO-50-1001 degrees C 

For a value of nu=17.06e-06 [13.6e-06 -> 23.06e-06] (m 2 /s), the uncertainty in the velocity prediction was 0 - 
27.727 percent as shown in Figure 7. 



0000 


Figure 7: Kinematic Viscosity - Velocity Uncertainty Percentage 


Gridsize 

For a grid size of 1,192,000 cells [grid 2 -1,862,500 cells], [grid3 - 3,31 1,689 cells], the uncertainty in the velocity 
prediction was 0 - 698 percent as shown in Figure 8. 



Figure 8: Grid Size - Velocity Uncertainty Percentage 

Turbulence Models 

The ke-realiable, kwSST, and SpalartAllmaras turbulence models converged using OpcnFoam and the uncertainty 
was calculated as an oscillatory input parameter as shown in Figure 9 



Figure 9: Turbulence Models - Velocity Uncertainty Percentage 

Solver 

OpenFoan, and Fluent were used to calculate the velocity distribution on the backward facing step and the 
uncertainty was calculated as an oscillatory input parameter as shown in Figure 10. 



Figure 10: Solver - Velocity Uncertainty Percentage 


The uncertainties of rhe variables with monotonic convergence (numerical) are calculated using Richardson's 
extrapolation as outlines by ASME V&V-20O9 «. This is accomplished Ihrough the five-step procedure Step I, 
calculate representative grid size, h as shown in equation 4. 


h . = 1 

f Total Volume \3 

1 1 
h^-( 

[total number of cells in fine grid) 
Total Volume 

\total number of cells in medium grid < 

ho = 1 

f Total Volume \? 

“3 | 

[total number of cells in coarse grid) 


( 4 ) 



Step 2 is to select three significantly (r>l .3) grid sizes and computer the ratio as shown i 


in equation 5. 


r 32 = 


A. 

_ *3 


( 5 ) 


Step 3 is to calculate the observed order, p, as shown in equation 6. This equation must be solved iteratively. 

P = l~\* [In fa) + In 

Un(r 21 )J \c 21 J \r 3 2 p -sign(^jj (®) 

Step 4 is to calculate the extrapolated values as shown in equation 7. 


c 21 _ ( r 2 i p - S k2 ) 

°ext — Z — 

0" 2 i p - 1) 


, 21 _ (S kl -S k2 ) 

a ~ (7> 


Step 5 is to calculate the fine grid convergence index and numerical uncertainty as shown in equation 8. This 

approached used a factor of safety of 1.25 and assumed that the distribution is Gaussian about the fine grid, 90 % 
confidence. 


GCI fine 21 = 


1.25 » e a 21 

('-2i p - 1) 


U 


mono tonic ~~ 


_ Mlfine 21 


1.65 


( 8 ) 


Numerical 


For a grid 
prediction 



size of 1,192,000 cells [grid 2 -1,862,500 cells], [grid3 - 3,31 1,689 cells], the uncertainty in the velocity 
was 0 - 5300 percent as shown in Figure 1 1 as estimated by Richardson’s extrapolation method. 


Figure 1 1 Numerical - Velocity Uncertainty Percentage 


A root-sum-squared (rss) of the uncertainty variables was calculated (omitting Richardson’s Extrapolation - see 
Discussion) and the velocity magnitude is shown in figure 12 with the corresponding uncertainty. 
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Figure 12: Velocity Prediction and Uncertainty Plot for ke-realizabie Turbulence Model 
Tlte highest uncertainty is +/- 4 85 m/s. This occurs in the region shown in Figure 13 in red. Figure 13 is the 
same data presented on the right hand side of Figure 12, except zoomed in to the region near the backward step and 


a smaller scale is used. 
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Figure 13: Velocity Uncertainty Plot for ke-realizable Turbulence Model 


V. Discussion 


The monotonic convergence uncertainty calculation (equations 4-8) was omitted in the rss uncertainty plot due to 
the values that were produced by using this method. The method produced uncertainty values that were on the order 
of 5000 percent of the localized velocity in the region near the backward step. It is believed this is due to the 
turbulence. Turbulence is calculated as a steady state value and fluctuations about that steady state. The 
fluctuations are inducing a non-linear result between the three grids and providing very large uncertainty bands in 
the localized region near the backward step. However, once you move approximately 5 lengths downstream of the 
backward step, the method begins producing reasonable results of 0 - 30 percent of the localized velocity. Treating 
the highly turbulent region behind the backward step as a monotonic case is in appropriate. It is believed that 
treating the grid as an input parameter with oscillatory convergence provides better results for a steady state, 
turbulent CFD simulation. This is evident in the values shown in Figure 1 . Most of the cells are exhibiting 
oscillatory convergence. It is believed all cells are exhibiting oscillatory convergence, however depending on when 
the sample takes place, one could misrule the results as monotonic or divergent. 

During several preliminary cases of the grid convergence study, one case provided an excellent example of 
domain sizing. A CFD analyst is always troubled with trying to keep the domain size large enough to not affect the 
solution. Using the oscillatory method, one can see the solution differences between the three grids. In the case 
presented below, the domain size is too small. This is evident by calculating the uncertainty using an oscillatory 
convergence method as shown in Figure 14. 



Figure 14: Example of (3) Grids with a Domain not large enough for the calculation 


The uncertainty is high near the boundary conditions; this shows the domain is not sized appropriately. 

It is computationally time consuming to run all of the different input parameters. It is suggested that the 
community compile a list of input parameters for each of the turbulence models and estimated uncertainty values for 
each of the parameters. The code used for this study was OpenFOAM and Table 2 of the Appendix is a list of all the 
input parameters for a ke-realizable case. Table 2 includes the input parameters and estimated uncertainties for the 
backward facing step. All values presented are in the form of a percent of the localized velocity. 

VI. Conclusion 

This paper outlines an uncertainty analysis for the ke realizable turbulence model for a backward facing step. 
The velocity magnitude was predicted using CFD. The uncertainty parameters listed in Table 1 were analyzed using 
an oscillatory convergence calculation or a monotonic convergence calculation. Plots of the velocity magnitude can 
be combined with a corresponding uncertainty plot for an accurate velocity prediction. 

There are other variables that would influence the uncertainty calculation. Examples of these other parameters 
include solution schemes, other turbulence models, and time accurate solutions. Future work will include analyzing 
each of these items. 

It is suggested that the CFD community begin to compile a list of the many variables associated with each 
uncertainty calculation for different problems and output variables. Ideally, an analyst could assemble a table of all 
uncertainty variables and estimate a number based on historical data rather than running separate CFD cases for each 
variable. The procedure above has been scripted and future work will include other geometries and turbulence 
models. 
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Table 2: Table of Input Parameters for ke-realizable OpenFOAM Turbulence Model 


ke Inputs: 





Default Values 

Recommended Uncertainty 








Boundary Conditions: 








epsilon 

inlet 

turbulentMixinglengthOissipationRatelnlet 

mixing length 

0.5 

1.20% 





value 

1 




wall 

epsilonWallFunction 

value 

0 




outlet 

type 

value 

1 










k 

inlet 

turbulentlntensityKineticEnergYlnlet 

intensity 

0.05 

0.80% 





value 

1 











wall 

kqRWallFunction 

value 

0 











outlet 

type 

value 

1 










nut 

inlet 

type 

value 

1 




wall 

nutWallFunction 

value 

0 




outlet 

type 

uniform 

0 










nuTilda 

inlet 

zeroGradient 






wall 

zeroGradient 






outlet 

zeroOradient 












P 

inlet 

zeroGradient 



lOx the variation 



wall 

zeroGradient 






outlet 

type 

fixed value 

0 










U 

inlet 

type 

X 

10 

1.3x the variation 





y 

0 






z 

0 











wall 

type 

X 

0 






y 

0 






z 

0 











outlet 

type 

X 

1 






y 

0 






z 

0 









Turbulence 

Turbulence Tranport Model 

nu 

Newtonian 

(02 100001 

1.00E-06 

(0-100 deg C) -> 28% 


Turbulence Propoerties 

simulationType 

RASModel 





RAS Properties 


RASModel 

realizableKE 










ControlDict 

End Time 




150000 



timestep (or CFL) 




1 



write percision 




6 



timePrecision 




6l 



Table 2: Table of Input Parameters for ke-realizable OpenFOAM Turbulence Model - continued 


SolutionSchemes 

gradSchemes 

P 

Gauss linear 






U 

Gauss Linear 












divSchemes 

phi, U 

Gauss limitedbnearV 






phi,k 

Gauss limitedlinear 






phi, epsilon 

Gauss limitedlinear 






phi,R 

Gauss limitedlinear 






R 

Gauss linear 






phi, nuTilda 

Gauss limitedlinear 






nuEff*dev(T(grad(U))) 

Gauss linear 












laplacianSchemes 

nuEff,U 

Gauss linear corrected 






1| A(U)),p 

Gauss linear corrected 






DkEff.k 

Gauss linear corrected 






DepsilonEff, epsilon 

Gauss linear corrected 






DREff.R 

Gauss linear corrected 






DnuTildaEff, nuTilda 

Gauss linear corrected 












interpolationSchemes 

interpolate(U) 

linear 












snGradSchemes 

default 

corrected 





fluxRequired 

default no 

P 











Solvers 

P 

solver 

GAMG; 



30% 



smoother 

GaussSeidel; 






cache Agglomeration 

true; 






nCellstnCoarsestLevel 

10; 






agglomerator 

faceAreaPair; 






mergeLevels 

1; 






tolerance 

le-06; 






relTol 

0.05; 












pFinal 

solver 

GAMG; 






smoother 

GaussSeidel; 






cacheAgglomeration 

true; 






nCellsInCoarsestLevel 

10; 






agglomerator 

faceAreaPair; 






mergeLevels 

1; 






tolerance 

le-06; 






relTol 

0; 












(U|k|epsilon)Fmal 

solver 

PBiCG; 






preconditioner 

DILU; 






tolerance 

le-05, 






relTol 

0; 












PIMPLE 

nOuterCorrectors 

4; 






nCorrectors 

1; 






nNonOrthogonalCorrectors 

0; 






pRefCell 

0; 






pRefValue 

0; 












SIMPLE 

nNonOrthogonalCorrectors 

0; 






residualControl 

P 

le-2; 






U 

le-3; 






(k | epsilon) 

le-3; 











relaxationFactors 

P 

0.3; 






U 

0.7; 






k 

0.7; 






epsilon.* 

0.7; 





cache 

grad(U); 












Mesh Size 

number of cells 





Case Specific - Use Oscillatory Method 
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Background 


• How good is CFD? 



• There are uncertainties and errors in using CFD 

- No standard method for evaluating uncertainty in CFD 

- Potential Errors include: 


• physical approximation error 

• computer round-off error 

• iterative convergence error 

• discretization errors 

• computer programming errors 

• usage errors 

• turbulence induced errors 
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Background 




• This presentation describes an approach to 
validate the uncertainty in using CFD. The 
method will use the state of the art 
uncertainty analysis applied to the ke- 
realizable turbulence model to predict the 

velocity uncertainty of a backward facing 
step. 
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Velocity Magnitude Prediction - Backward 

Facing Step 



LAUNCH SERVICES 


Symmetry 
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Velocity Prediction with Uncertainty 



* 

#>UCF J 


Uncertainly in Velocity (Uvel) ♦/- (m/sj 

■ 2 45 
2 30 
2.14 
1.99 
1 84 
1 68 
1 53 
138 

1.22 Wf*'* 

1.07 
0 92 
t 0 77 
0 61 
> 0 46 
0 31 
0 15 
0 00 
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Method 



• To estimate the uncertainty, the following 
ASME Standard was used. 


ASME V&V 20-2009 "Standard for Verification and 
Validation in Computational Fluid Dynamics and 

Heat Transfer" 
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Summary of Method 



• Convergence studies require a minimum of three solutions to 
evaluate convergence with respect to an input parameter. 
Consider the situation for 3 solutions corresponding to fine S kl , 
medium S k2 , and coarse S k3 values for the kth input parameter. 
Solution changes e for medium-fine and coarse-medium solutions 
and their ratio R k are defined by: 

£ 21 = ^k2 ~~ ^kl 
£ 32 = ^k3 ~~ ^k2 
^k = £ 21 / E 32 

• Three convergence conditions are possible: 

Monotonic convergence: 0< R k <1 
Oscillatory convergence: R k < O' 

Divergence: R k >l 
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• The uncertainty associated with the CFD calculation is the compilation of 
the elemental errors associated with each of the numerical, input, and 
solver errors. This uncertainty can be calculated using a Data Reduction 
equation the form r = r(X v X 2 ,...Xj) as shown, 


= (el {(f ) 2 Bf] + 2 z,LzLi+i {(f) (f) r B,B k ] correlated ) + ZL j(f ) 2 /»,=*}) 


l / 2 


• Where, 

B| 

(®i ^k)correlated “ 

and Pj 


the systematic (bias) error associated with variable Xj, 

the correlated systematic error between variables Xjand X k , 
the random error associated with variable Xj. 


For the calculation, the correlated errors and random errors are neglected 
and the data reduction equation reduces to the following, as shown 
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Uncertainty Variables ke-realizable 
(OPENFOAM - SimpleFoam) 



LAUNCH SERVICES 


• There are 87 Different Input Parameters for the 
ke-realizable model in SimpleFoam 
— These include: 


• Boundary Conditions 

• Wall Functions 

• Fluid Properties 

• Turbulence Parameters 

• Solution Schemes 

• Solvers 

• Mesh 

• ect. 
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Uncertainty Variables Considered 
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Type of Variable 

Variables Xi 

Value 

Bias Error 

Boundary 

Conditions 

epsilion turbulent mixing length dissipation rate inlet (m2/s3) 

0.5 

0.5 


k turbulent intensity kinetic energy inlet (m2/s2) 

0.05 

0.05 


pressure outlet (Pa) 

101325 

2 % 


velocity inlet (m/s) 

10 

0.5 

Fluid Properties 

kinematic viscosity nu represents air (0-50-100) deg C 

1.79E-06 

[13.6e-06-> 23.06e- 
06] 

Grid Size 

Method - Uses Oscillatory Uncertainty 

1,192.000 

1,862,500 

3,311,689 


Numerical 

Method - Uses Richardson’s Extrapolation (ASME 5 Step Procedure) - Calculated for Velocity at each 

Cell 

Solver 

OpenFOAM (SimpleFoam) vs. Fluent 



Turbulence 

Models 

ke-realiable, kwSST, and SpalartAllmaras 
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Expanding the data reduction equation for the listed variables in order from top 
to bottom. 




Oscillatory Variables 
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• The uncertainty for each of the following was calculated for each 
cell using the following method outlined by Stern, Wilson, Coleman, 
and Paterson. S is the simulated result. For this case it is the upper 
velocity S u and the lower velocity S L . 

- epsilion turbulent mixing length dissipation rate inlet (m 2 /s 3 ) 

- k turbulent intensity kinetic energy inlet (m 2 /s 2 ) 

- Pressure outlet (Pa) 

- Velocity Inlet (m/s) 

- Kinematic viscosity nu=17.06e-06 f!3.6e-06 -> 23.06e-061 (m 2 / 
s) represents air [0-50-1001 degrees C 

- Grid size 

— Turbulence Models 

— Solver 



U Oscillatory 2 ^l) 
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Results (Oscillatory Variables) 
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LAUNCH SERVICES PROGRAM 


epsilion turbulent mixing 
length diss ipation rate 
inlet (m 2 /s 3 ) 



Pressure outlet (Pal 

0-20 percent 



Kinematic viscosity 
nu=17.06e-06 fl3.6e-06 -> 
23.06e-061 (m 2 /s) represents 
air [0-50-10 01 degrees C 

0-27.727 percent 



Lx 

Turbulence Models 

> 100 % 



k turbulent intensity 
kinetic energy inlet (m 2 / 
sfl 

0-0.785 percent 



Velocity Inlet (m/si 

0-6.558 percent 




Grid size 

0 - 698 percent 


Solver 

> 30 % 




Percent - is the percentage change in local velocity 
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Monotonic Convergence Variables 

(Numerical) 

D ' • 
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• The uncertainties of the variables with monotonic 
convergence (numerical) are calculated using Richardson's 
extrapolation as outlines by ASME V&V-2009. This is 
accomplished through the five-step procedure. 

• Step 1, calculate representative grid size, h as shown 



total number of cells in fine grid 



total number of cells in medium grid 


i 
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Monotonic Convergence Variables 
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Step 2, calculate representative grid ratio, r as shown 


r 21 — 


K 



• Step 3 is to calculate the observed order, p, as shown. This 
equation must be solved iteratively. 
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Monotonic Convergence Variables 

(Numerical) 


LAUNCH 



• Step 4 is to calculate the extrapolated values as shown 

c 21 ( r 2l P * 1 — S fc2 ) 


'ext ~ 


(r 2 i p ~ 1) 


p 21 ($ki £k2} 

0 ~ (Ski) 

Step 5 is to calculate the tine grid convergence index and 
numerical uncertainty as shown. This approached used a 
factor of safety of 1.25 and assumed that the distribution is 
Gaussian about the fine grid, 90 % confidence. 
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GCl fi 


21 


ne 


1.25 » e a 21 
(r 2 iP - 1) 


U 


monotonic 


GCI 


fine 


1.65 


21 


h 


Results (Monotonic 
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Numerical 



For a grid size of 1,192,000 cells [grid 2 -1,862,500 
cells], [grid3 - 3,311,689 cells], the uncertainty in 
the velocity prediction was 0 - 5300 percent as 

shown in Figure 11 as estimated by Richardson's 
extrapolation method. 
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Results (Monotonic 
Convergence) 

Numerical 

Three convergence conditions are possible: 

Monotonic convergence: 0< R k <1 
Oscillatory convergence: R k < O' 
Divergence: R k >l 



It is believed the errors in this method is due to the turbulence. Turbulence is calculated as a 


steady state value and fluctuations about that steady state. The fluctuations are inducing a non- 
linear result between the three grids and providing very large uncertainty bands in the localized 
region near the backward step. Treating the highly turbulent region behind the backward step as 
a monotonic case is in appropriate. 
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A root-sum-squared (rss) of the uncertainty 
variables was calculated (omitting 
Richardson's Extrapolation) 



Uncertainity in Velocity (Uvel) ♦/- [m/sl 

H 4 85 
4 54 1 
4 24 
3 94 
363 
333 
r 3.03 
2.73 
2 42 
2 12 
1 82 
1 51 
1.21 
0 91 
061 
030 
0 00 
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The highest uncertainty is +/- 4.85 m/s. 



• During several preliminary cases of the grid 

convergence study, one case provided an excellent 
example of domain sizing. 

— A CFD analyst is always troubled with trying to keep the 
domain size large enough to not affect the solution. 

• Using the oscillatory method, one can see the solution differences 
between the three grids. In the case presented below, the domain 
size is too small. 
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Conclusion 

launch services program ■ 

• This paper outlines an uncertainty analysis for the ke ^ ' 
realizable turbulence model for a backward facing step. 

• The velocity magnitude was predicted using CFD. 

• The uncertainty parameters listed in Tablel were analyzed 
using an oscillatory convergence calculation or a monotonic 
convergence calculation. 

• Plots of the velocity magnitude can be combined with a 
corresponding uncertainty plot for an accurate velocity 
prediction. 

• There are other variables that would influence the 
uncertainty calculation. Examples of these other 
parameters include solution schemes, other turbulence 
models, and time accurate solutions. Future work will 
include analyzing each of these items. 
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Conclusion /Recommendation 
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LAUNCH 


The following input uncertainty's are recommended 



Type of 
Variable 

Variables Xi 

Value 

Bias Error 

Uncertainly 

Boundary 

Conditions 

epsilion turbulent mixing length dissipation rate inlet (m2/s3) 

0.5 

0.5 

1.2% of local velocity 


k turbulent intensity kinetic energy inlet (m2/s2) 

0.05 

0.05 

0.8 % of local velocity 


pressure outlet (Pa) 

101325 

2% 

lOx the variation 


velocity inlet (m/s) 

10 

0.5 

1.3x the variation 

Fluid 

Properties 

kinematic viscosity nu represents air [0-50-100] deg C 

1.79E-06 

[13.6e-06 -> 
23.06e-06] 

28% of the local velocity 

Grid Size 

Method - Uses Oscillatory Uncertainty 

1,192,000 

1,862,500 


grid specific 

3,311,689 


Numerical 

Method - Uses Richardson's Extrapolation (ASME 5 Step 
Procedure) - Calculated for Velocity at each Cell 



Solver 

OpenFOAM (SimpleFoam) vs. Fluent 



30% of the local velocity 

Turbulence 

Models 

ke-realiable, kwSST, and SpalartAllmaras 



Future work will 
consider more 
turbulence models 



UCF 


University of Central Florida 


COLLEGE OF ENGINEERING AND COMPUTER SCIENCE 


JOHN f KENNEDY SPACE CENTER 


Recommendation / Future Work 



• It is suggested that the CFD community begin 
to compile a list of the many input parameters 
associated with each uncertainty calculation 
for different problems and output variables. 

• Ideally, an analyst could assemble a table of 
all uncertainty variables and estimate a 
number based on historical data rather than 
running separate CFD cases for each variable. 
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Backup 
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Uncertainty Variables ke-realizable 
(OPENFOAM - SimpleFoam) 



LAUNCH SERVICES PROGRAM 


ke Inputs: 





Default Values 







Boundary Conditions: 







epsilon 

inlet 

turbulent Mixing Length Dissipation Rate In let 

mixing length 

0.5 





value 

1 



wall 

ep si lonWa II Function 

value 

O 



outlet 

type 

value 

1 








k 

inlet 

turbulentlntensityKineticEnergylnlet 

Intensity 

0.05 





value 

1 









wall 

kqR Wall Function 

value 

o 









outlet 

type 

value 

1 








nut 

inlet 

type 

value 

1 



wall 

nutWallFunction 

value 

o 



outlet 

type 

uniform 

o 








nuTilda 

inlet 

zeroGradient 





wall 

zeroG radient 





outlet 

zeroGradient 










P 

inlet 

zeroGradient 





wall 

zeroGradient 





outlet 

type 

fixed value 

o 








U 

inlet 

type 

X 

io 





V 

o 





z 

o 









wall 

type 

X 

o 





V 

o 





z 

o 









outlet 

type 

X 

1 





V 

o 





z 

o 







Turbulence 

Turbulence Tranport Model 

nu 

Newtonian 

(02-10000] 

1.00E-06 


Turbulence Propoerties 

simulationType 

RASModel 




RAS Properties 


RASModel 

realizableKE 








ControlDict 

End Time 




150000 


timestep (or CFL) 




1 


write percision 




6 


timePrecision 




6 
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SolutlonSchemes 

gradSchemes 

P 

Gauss linear 




U 

Gauss Linear 








divSchemes 

phi, LJ 

Gauss limitedLinearV 




phi.k 

Gauss IlmitedLinear 




phl.epsllon 

Gauss IlmitedLinear 




phi, R 

Gauss IlmitedLinear 




R 

Gauss linear 




phi, nuTilda 

Gauss IlmitedLinear 




nuEff*dev(T(grad(U))) 

Gauss linear 








laplaclanSchemes 

nu Eff , LJ 

Gauss linear corrected 




l|A(U)),p 

Gauss linear corrected 




DkEff, k 

Gauss linear corrected 




Deps i Ion Eff , epsilon 

Gauss linear corrected 




DREff, R 

Gauss linear corrected 




DnuTlldaEff,nuTllda 

Gauss linear corrected 








InterpolationSchemes 

lnterpolate(U ) 

linear 








snGradSchemes 

default 

corrected 



fluxRequired 

default no 

P 







Solvers 

£ 

solver 

GAMG; 




smoother 

GaussSeidel; 




cache Agglomeration 

true; 




nCel Is InCoarsest Level 

lO; 




«*8glomerator 

faceAreaPalr; 




merge Levels 

X; 




tolerance 

le-06; 




relTol 

0.05; 








pFInal 

solver 

GAM G ; 




smoother 

GaussSeidel; 




cache Agglomeration 

true; 




nCellsInCoarsest Leve 1 

XO; 




agglomerator 

faceAreaPalr; 




merge Leve Is 

X; 




tolerance 

Xe-06; 




relTol 

O; 








( LJ | k | epsilon) Final 

solver 

PBICG; 




preconditioner 

Dl LLJ ; 




tolerance 

Xe-OS; 




relT ol 

O; 








PIIVIPLE 

nOuterCorrectors 

■4; 




nCorrectors 

X; 




oNonOrthogonalCorrectors 

O; 




pRefCel 1 

O; 




pRefValue 

O; 








SIMPLE 

nNonOrthogonalCorrectors 

O; 




residualControl 

P 

Xe-2; 




U 

Xe-3; 




(k | epsilon) 

Xe-3 ; 







relaxation Factors 

£ 

0.3; 




LJ 

0 . 7 ; 




k 

0 . 7 ; 




epsilon. * 

0 . 7 ; 



cache 

grad(U); 








Mesh Size 

number of cells 
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